#!/bin/bash

source env.sh
source parameters.sh
source functions.sh
source $calj_init_conf/parameters.sh

cwd=`pwd`
mylog=$cwd/cal.j.log
rm -f $mylog

calj_epsilon=$calj_f
calj_sigma=`echo "sqrt( 2. / $calj_gamma )" | bc -l`

rm -fr tools
cp -a $calj_tool_dir ./tools
calj_tool_dir=$cwd/tools
make -C $calj_tool_dir/dihedral.table clean > /dev/null
make -j -C $calj_tool_dir/dihedral.table > /dev/null
make -C $calj_tool_dir/fmpt clean > /dev/null
make -j -C $calj_tool_dir/fmpt > /dev/null
calj_meta_width=`echo "180 - $calj_meta_low" | bc -l`
echo "# calj_meta_width is $calj_meta_width"

# make init table
$calj_tool_dir/dihedral.table/gmx.funct1 --kphi $calj_init_dih_k --phi $calj_init_dih_phi --multi $calj_init_dih_n --bin $calj_tab_bin --output table.v.xvg &>> $mylog
if test ! -f $calj_init_j_tab; then
    echo "# no init table, generate one"
    $calj_tool_dir/dihedral.table/gmx.funct1 --kphi $calj_init_j_k --phi $calj_init_j_phi --multi 1 --output table.j.init.xvg &>> $mylog
else
    echo "# use init table $calj_init_j_tab"
    mv $calj_init_j_tab ./table.j.init.xvg.tmptmptmp
    mv ./table.j.init.xvg.tmptmptmp ./table.j.init.xvg
fi
ppi=`printf %03d 0`
dir_name=step.$ppi
if test -d $dir_name; then
    mv $dir_name $dir_name.`date +%s`
fi
mkdir -p $dir_name
mv table.v.xvg  $dir_name
cp table.j.init.xvg  $dir_name/table.j.new.xvg
jj=0
# sqrt2gamma=`echo "sqrt(2. * $calj_gamma)" | bc -l`

# do loop
for ii in `seq 1 $calj_num_iter`;
do
    ppi=`printf %03d $ii`
    dir_name=step.$ppi
    ppj=`printf %03d $jj`
    last_dir_name=step.$ppj
    echo "## doing step $ppi"
    
    if test -d $dir_name; then
	mv $dir_name $dir_name.`date +%s`
    fi
    mkdir -p $dir_name
    cd $dir_name

    cp ../$last_dir_name/table.v.xvg .
    cp ../$last_dir_name/table.j.new.xvg ./table.j.init.xvg
    $calj_tool_dir/dihedral.table/combine.table --input-0 table.v.xvg --input-1 table.j.init.xvg --scale-0 1.0 --scale-1 2.0 --output table_d0.xvg
    
    cp $calj_gmx_template/grompp.mdp .
    cp $calj_gmx_template/topol.top .
    cp $calj_gmx_template/angle.ndx .
    test -f $calj_gmx_template/index.ndx && cp $calj_gmx_template/index.ndx .

    # cal trajs from init files
    rm -f trajs
    mkdir -p trajs
    rm -f angaver.name traj.name
    cd trajs
    for angle in `seq $init_interval_start $init_center_angle_interval $init_interval_end`
    do
	pa=`printf %+04d $angle`
	init_nstep=`echo "$init_prod_time / $init_frame_feq" | bc -l | cut -f 2 -d '=' | cut -f 1 -d '.'`
	echo "# doing angle $angle, $init_nstep initial confs"
	for step in `seq 1 $init_nstep`
	do
	    pstep=`printf %06d $step`
	    work_dir_name=traj.$pa.$pstep
	    rm -f $work_dir_name
	    mkdir -p $work_dir_name
	    cd $work_dir_name
	    ln -s ../../topol.top .
	    ln -s ../../angle.ndx .
	    test -f ../../index.ndx && ln -s ../../index.ndx .
	    ln -s ../../table_d0.xvg .
	    cp ../../grompp.mdp .
	    set_parameters grompp.mdp
	    my_time=`echo "($step - 0.5) * $init_frame_feq" | bc -l`
	    if test -f index.ndx; then
		echo "# grompp command: grompp -n index.ndx -c $calj_init_conf/init.conf.$pa/conf.gro -t $calj_init_conf/init.conf.$pa/traj.trr -time $my_time" &>> run.log
		grompp -n index.ndx -c $calj_init_conf/init.conf.$pa/conf.gro -t $calj_init_conf/init.conf.$pa/traj.trr -time $my_time &>> run.log
	    else
		echo "# grompp command: grompp -c $calj_init_conf/init.conf.$pa/conf.gro -t $calj_init_conf/init.conf.$pa/traj.trr -time $my_time" &>> run.log
		grompp -c $calj_init_conf/init.conf.$pa/conf.gro -t $calj_init_conf/init.conf.$pa/traj.trr -time $my_time &>> run.log
	    fi
	    mdrun  &>> run.log
	    if test $? -ne 0; then
		echo "# failed mdrun, exit"
		exit
	    fi
	    g_angle -type dihedral -od angdist.xvg  -ov angaver.xvg -xvg none &>> run.log
	    echo "./trajs/$work_dir_name/angaver.xvg" >> ../../angaver.name
	    echo "./trajs/$work_dir_name/traj.xtc"    >> ../../traj.name
	    rm -f state*cpt
	    if test $step -eq 1; then
		cd ..
		cp -a $work_dir_name backup.$work_dir_name
		cd $work_dir_name
	    fi
	    rm -f angdist.xvg angle.ndx ener.edr grompp.mdp *log table_d0.xvg topol.top mdout.mdp  topol.tpr confout.gro
	    cd ..
	done
    done
    cd ..
    echo "# cal j using command: $calj_tool_dir/fmpt/cal.j.recycle --meta-center 180 --meta-width $calj_meta_width --input-angle angaver.name --input-traj traj.name --input-j table.j.init.xvg --sigma $calj_sigma --epsilon $calj_epsilon --bin $calj_stat_bin --output fmpt.out --angle-dx $calj_grad_angle_dx"
    $calj_tool_dir/fmpt/cal.j.recycle --meta-center 180 --meta-width $calj_meta_width --input-angle angaver.name --input-traj traj.name --input-j table.j.init.xvg --sigma $calj_sigma --epsilon $calj_epsilon --bin $calj_stat_bin --output fmpt.out --angle-dx $calj_grad_angle_dx
    csg_resample --in fmpt.out --out v.out --derivative d.out --grid -180:$calj_tab_bin:180 --boundaries periodic
    $calj_tool_dir/dihedral.table/join.table -v v.out -d d.out --output table.j.new.xvg
    rm -f v.out d.out

    rm -fr trajs/traj.*
    cd $cwd
    jj=$ii
done

